#include "ProbabilityDistribution.h"

#include <gsl/gsl_math.h>
#include <gsl/gsl_errno.h>

#include <cstdio>
#include <iostream>
#include <limits>
#include <utility>
#include <sstream>
#include <boost\progress.hpp>

ProbabilityDistribution::ProbabilityDistribution(void) : _zeroQuantileLocation(-1.0), _mean(0.0), _variance(0.0)
{
	// initialize the solver
	_solver = gsl_root_fsolver_alloc (gsl_root_fsolver_brent);

	// initialize the integrator
	_integrationWorkspace = gsl_integration_workspace_alloc (1000);

}


ProbabilityDistribution::~ProbabilityDistribution(void)
{
	gsl_root_fsolver_free (_solver);
	gsl_integration_workspace_free (_integrationWorkspace);
}


double ProbabilityDistribution::cdf(double x) 
{
  
    double result, error;      
	gsl_function F;
	F.function = &ProbabilityDistribution::densityWrapper;
	F.params = this;
	
	double epsabs = 0;
	double epsrel = 1e-13;

	try{
		gsl_integration_qagil (&F, x, epsabs, epsrel, 1000, _integrationWorkspace, &result, &error);
	}catch(...)
	{
		std::stringstream errMsg;
		errMsg << "failed to evaluate cdf at x = " << x;
		throw std::runtime_error(errMsg.str());
	}
	//printf("result          = % .18f\n", result);
    //printf ("estimated error = % .18f\n", error);
    //printf ("intervals =  %d\n", w->size);
     
	return result;

}

double ProbabilityDistribution::quantile(double u) 
{
	//std::cout << "quantile evaluation at u = " << u << std::endl;
	//boost::progress_timer t;

	int status;
	double epsrel = 1e-15;
    int iter = 0, max_iter = 1000;
        
    double r = 0; // root

	std::pair<double, double> qBounds = quantileBounds(u);
    double x_lo = qBounds.first, x_hi = qBounds.second; // interval enclosing root
    
	RootInfo info;
	info.distribution = this;
	info.u = u;

	gsl_function F;     
	F.function = &ProbabilityDistribution::cdfWrapper;
    F.params = &info;
     
    gsl_root_fsolver_set (_solver, &F, x_lo, x_hi);
     
    //printf ("using %s method\n", gsl_root_fsolver_name (s));
    //printf ("%5s [%9s, %9s] %9s %10s %9s\n", "iter", "lower", "upper", "root", "err", "err(est)");
    
	try{
		do{
			iter++;
			status = gsl_root_fsolver_iterate (_solver);
			r = gsl_root_fsolver_root (_solver);
			x_lo = gsl_root_fsolver_x_lower (_solver);
			x_hi = gsl_root_fsolver_x_upper (_solver);
			status = gsl_root_test_interval (x_lo, x_hi, 0, epsrel);
     
			//if (status == GSL_SUCCESS)
			//    printf ("Converged:\n");
     
			//printf ("%5d [%.7f, %.7f] %.10f %+.7f %.7f\n", iter, x_lo, x_hi, r, x_hi - x_lo);
		}
		while (status == GSL_CONTINUE && iter < max_iter);
	}catch(...){
		std::stringstream errMsg;
		errMsg << "quantile evaluation: root finder failed at u = " << u;
		throw std::runtime_error(errMsg.str());
	}

    return r;
}

double ProbabilityDistribution::zeroQuantileLocation()
{
	_zeroQuantileLocation = (_zeroQuantileLocation == -1) ? cdf(0.0) : _zeroQuantileLocation;
	return _zeroQuantileLocation;
}

